home *** CD-ROM | disk | FTP | other *** search
/ The Arsenal Files 4 / The Arsenal Files 4 (Arsenal Computer).ISO / ham / sattrk31.tgz / sattrack-3.1.tar / SatTrack / src / sattrack / satprop.c < prev    next >
C/C++ Source or Header  |  1995-03-16  |  27KB  |  739 lines
  1. /******************************************************************************/
  2. /*                                                                            */
  3. /*  Title       : satprop.c                                                   */
  4. /*  Author      : Manfred Bester                                              */
  5. /*  Date        : 20Nov93                                                     */
  6. /*  Last change : 15Mar95                                                     */
  7. /*                                                                            */
  8. /*  Synopsis    : SGP4/SDP4 propagation routines (translated from FORTRAN)    */
  9. /*                for reference see:                                          */
  10. /*                "Spacetrack Report No. 3, Models for Propagation of         */
  11. /*                NORAD Element Sets", F. R. Hoots and R. L. Roehrich,        */
  12. /*                December 1980                                               */
  13. /*                                                                            */
  14. /*  REMARK: The SDP4 model implementation has not been completed in this      */
  15. /*          release of SatTrack. The full implementation of the SDP4 model    */
  16. /*          for deep space objects is expected to be part of a future         */
  17. /*          release. For now, a simplified model is used for deep space       */
  18. /*          objects (the so-called 'TLE Mean model').                         */
  19. /*                                                                            */
  20. /*                                                                            */
  21. /*  SatTrack is Copyright (c) 1992, 1993, 1994, 1995 by Manfred Bester.       */
  22. /*  All Rights Reserved.                                                      */
  23. /*                                                                            */
  24. /*  Permission to use, copy, and distribute SatTrack and its documentation    */
  25. /*  in its entirety for educational, research and non-profit purposes,        */
  26. /*  without fee, and without a written agreement is hereby granted, provided  */
  27. /*  that the above copyright notice and the following three paragraphs appear */
  28. /*  in all copies. SatTrack may be modified for personal purposes, but        */
  29. /*  modified versions may NOT be distributed without prior consent of the     */
  30. /*  author.                                                                   */
  31. /*                                                                            */
  32. /*  Permission to incorporate this software into commercial products may be   */
  33. /*  obtained from the author, Dr. Manfred Bester, 1636 M. L. King Jr. Way,    */
  34. /*  Berkeley, CA 94709, USA. Note that distributing SatTrack 'bundled' in     */
  35. /*  with ANY product is considered to be a 'commercial purpose'.              */
  36. /*                                                                            */
  37. /*  IN NO EVENT SHALL THE AUTHOR BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT, */
  38. /*  SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OF   */
  39. /*  THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE AUTHOR HAS BEEN ADVISED  */
  40. /*  OF THE POSSIBILITY OF SUCH DAMAGE.                                        */
  41. /*                                                                            */
  42. /*  THE AUTHOR SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT      */
  43. /*  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A   */
  44. /*  PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"      */
  45. /*  BASIS, AND THE AUTHOR HAS NO OBLIGATIONS TO PROVIDE MAINTENANCE, SUPPORT, */
  46. /*  UPDATES, ENHANCEMENTS, OR MODIFICATIONS.                                  */
  47. /*                                                                            */
  48. /******************************************************************************/
  49.  
  50. #include <stdio.h>
  51. #include <math.h>
  52.  
  53. #ifndef STDLIB
  54. #include <stdlib.h>
  55. #endif
  56.  
  57. #include "satglobalsx.h"
  58. #include "satglobalspx.h"
  59. #include "sattrack.h"
  60. #include "satprop.h"
  61.  
  62. /******************************************************************************/
  63. /*                                                                            */
  64. /* checkPropModelType: checks type of propagation model used                  */
  65. /*                                                                            */
  66. /******************************************************************************/
  67.  
  68. void checkPropModelType()
  69.  
  70. {
  71.     getSemiMajorAxis();
  72.     checkPerigeeHeight();
  73.  
  74.     if (deepSpaceFlag)
  75.     {
  76.         sprintf(propModel,"%s",TLEMEANMODEL);                 /* for now only */
  77.         propModelType = TLEMEAN;
  78.  
  79. /*
  80.         sprintf(propModel,"%s",DEEPSPACEMODEL);
  81.         propModelType = DEEPSPACE;
  82. */
  83.     }
  84.  
  85.     else
  86.     {
  87.         sprintf(propModel,"%s",LOWEARTHMODEL);
  88.         propModelType = LOWEARTH;
  89.     }
  90.  
  91.     /* 'modelFlag' is TRUE if main() is invoked with the '-m' option */
  92.     /* this overwrites the usage of either the SGP4 or SDP4 models   */
  93.  
  94.     if (modelFlag)  
  95.     {
  96.         sprintf(propModel,"%s",TLEMEANMODEL);
  97.         propModelType = TLEMEAN;
  98.     }
  99.  
  100.     return;
  101. }
  102.  
  103. /******************************************************************************/
  104. /*                                                                            */
  105. /* changeToNoradUnits: changes units of fundamental elements for use with     */
  106. /*                     NORAD model                                            */
  107. /*                                                                            */
  108. /******************************************************************************/
  109.  
  110. void changeToNoradUnits()
  111.  
  112. {
  113.     epoch0     = epochDay;                                   /* [d]           */
  114.  
  115.     ecc0       = eccentricity;                               /* [---]         */
  116.     inc0       = inclination;                                /* [rad]         */
  117.     meanMot0   = epochMeanMotion * TWOPI / MPD;              /* [rad/min]     */
  118.     meanAn0    = epochMeanAnomaly;                           /* [rad]         */
  119.     argPer0    = epochArgPerigee;                            /* [rad]         */
  120.     raan0      = epochRaan;                                  /* [rad]         */
  121.  
  122.     decRate    = decayRate * TWOPI / MPD2;                   /* [rad/min^2]   */
  123.     decRateDot = decayRateDot * TWOPI / MPD3;                /* [rad/min^3]   */
  124.  
  125.     bStar      = bStarCoeff / AE;
  126.  
  127.     return;
  128. }
  129.  
  130. /******************************************************************************/
  131. /*                                                                            */
  132. /* getSemiMajorAxis: recover original mean motion (meanMotDeep) and semimajor */
  133. /*                   axis (sma0Deep) from input elements                      */
  134. /*                                                                            */
  135. /******************************************************************************/
  136.  
  137. void getSemiMajorAxis()
  138.  
  139. {
  140.     changeToNoradUnits();
  141.  
  142.     cosInc        = cos(inc0);
  143.     sinInc        = sin(inc0);
  144.  
  145.     theta         = cosInc;
  146.     theta2        = SQR(theta);
  147.     theta4        = SQR(theta2);
  148.  
  149.     x1mth2        =  1.0 - theta2;
  150.     x3th2m1       =  3.0 * theta2 - 1.0;
  151.     x5th2m1       =  5.0 * theta2 - 1.0;
  152.     x7th2m1       =  7.0 * theta2 - 1.0;
  153.     x7th2m3       =  7.0 * theta2 - 3.0;
  154.     x19th2m4      = 19.0 * theta2 - 4.0;
  155.  
  156.     eccSQ         = SQR(ecc0);
  157.     eccCB         = ecc0 * eccSQ;
  158.  
  159.     beta02        = 1.0 - eccSQ;                          /* needed elsewhere */
  160.     beta0         = sqrt(beta02);                         /* needed elsewhere */
  161.     beta04        = SQR(beta02);                          /* needed elsewhere */
  162.  
  163.     sma1          = pow((KE/meanMot0),TWOTHIRDS);
  164.  
  165.     deltaFact     = 1.5 * CK2 * x3th2m1 / pow(beta02,THREEHALFS);
  166.     delta1        = deltaFact / SQR(sma1);
  167.     delta12       = SQR(delta1);
  168.     delta13       = delta1 * delta12;
  169.  
  170.     sma0          = sma1 * (1.0 - ONETHIRD * delta1 - delta12 
  171.                                 - 134.0/81.0 * delta13);
  172.  
  173.     delta0        = deltaFact / SQR(sma0);
  174.     meanMotDeep   = meanMot0 / (1.0 + delta0);
  175.     sma0Deep      = sma0 / (1.0 - delta0);
  176.     semiMajorAxis = sma0Deep * EARTHRADIUS;
  177.  
  178.     deepSpaceFlag = (TWOPI / meanMotDeep >= 225.0) ? TRUE : FALSE;
  179.  
  180.     return;
  181. }
  182.  
  183. /******************************************************************************/
  184. /*                                                                            */
  185. /* getTimeArgs: gets time arguments for the SGP4/SDP4 model calculations      */
  186. /*                                                                            */
  187. /******************************************************************************/
  188.  
  189. void getTimeArgs(tmArg)
  190.  
  191. double tmArg;
  192.  
  193. {
  194.     tFP = (tmArg - epoch0) * MPD;                                    /* [min] */
  195.     tSQ = tFP * tFP;
  196.     tCB = tSQ * tFP;
  197.     tQD = tCB * tFP;
  198.     tQN = tQD * tFP;
  199.  
  200.     return;
  201. }
  202.  
  203. /******************************************************************************/
  204. /*                                                                            */
  205. /* checkPerigeeHeight: checks perigee height and makes modifications          */
  206. /*                     accordingly                                            */
  207. /*                                                                            */
  208. /* For a perigee below 156 km, the values of S and Q0MS0T4 are altered.       */
  209. /*                                                                            */
  210. /******************************************************************************/
  211.  
  212. void checkPerigeeHeight()
  213.  
  214. {
  215.     double perigeeTest;
  216.  
  217.     perigeeTest   = (sma0Deep * (1.0 - ecc0)/AE - AE) * EARTHRADIUS;
  218.     perigeeHeight = (sma0Deep * (1.0 - ecc0)    - AE) * EARTHRADIUS;
  219.     apogeeHeight  = (sma0Deep * (1.0 + ecc0)    - AE) * EARTHRADIUS;
  220.  
  221.     truncFlag     = (perigeeTest   < 220.0) ? TRUE : FALSE;
  222.     satCrashFlag  = (perigeeHeight <   0.0) ? TRUE : FALSE;
  223.  
  224.     sStar = S;
  225.     qmst4 = Q0MS0T4;
  226.  
  227.     if (perigeeHeight < 156.0)
  228.     {
  229.         sStar = perigeeHeight - S0;
  230.  
  231.         if (perigeeHeight <= 98.0)
  232.             sStar = 20.0;
  233.     
  234.         qmst4 = pow(((Q0 - sStar) * AE / EARTHRADIUS),4.0);
  235.         sStar = sStar / EARTHRADIUS + AE;
  236.     }
  237.  
  238.     return;
  239. }
  240.  
  241. /******************************************************************************/
  242. /*                                                                            */
  243. /* calcCommonPerturb: calculates quantities that are common to the SGP4 and   */
  244. /*                    SDP4 models                                             */
  245. /*                                                                            */
  246. /******************************************************************************/
  247.  
  248. void calcCommonPerturb()
  249.  
  250. {
  251.     pInvSQ = 1.0 / (SQR(sma0Deep) * beta04);
  252.     xi     = 1.0 / (sma0Deep - sStar);
  253.     eta    = sma0Deep * ecc0 * xi;
  254.     etaSQ  = SQR(eta);
  255.     eccEta = ecc0 * eta;
  256.     psiSQ  = fabs(1.0 - etaSQ);
  257.     coef0  = qmst4 * pow(xi,4.0);
  258.     coef1  = coef0 / pow(psiSQ,3.5);
  259.  
  260.     c2     = coef1 * meanMotDeep * 
  261.              (sma0Deep * (1.0 + 1.5 * etaSQ + eccEta * (4.0 + etaSQ)) + 
  262.              0.75 * CK2 * xi / psiSQ * x3th2m1 * 
  263.              (8.0 + 3.0 * etaSQ * (8.0 + etaSQ)));
  264.  
  265.     c1     = bStar * c2;
  266.  
  267.     c3     = coef0 * xi * A3OCK2 * meanMotDeep * AE * sinInc / ecc0;
  268.  
  269.     c4     = 2.0 * meanMotDeep * coef1 * sma0Deep * beta02 * 
  270.              ((2.0 * eta * (1.0 + eccEta) + 0.5 * ecc0 + 0.5 * eta * etaSQ) - 
  271.              2.0 * CK2 * xi / (sma0Deep * psiSQ) * 
  272.              (-3.0 * x3th2m1 * (1.0 + 1.5 * etaSQ - 2.0 * eccEta - 
  273.              0.5 * eccEta * etaSQ) + 0.75 * x1mth2 * 
  274.              (2.0 * etaSQ - eccEta - eccEta * etaSQ) * cos(2.0 * argPer0)));
  275.  
  276.     c5     = 2.0 * coef1 * sma0Deep * beta02 * 
  277.              (1.0 + 2.75 * (etaSQ + eccEta) + eccEta * etaSQ);
  278.  
  279.     clc1   = 3.0  * CK2 * pInvSQ * meanMotDeep;
  280.     clc2   = clc1 * CK2 * pInvSQ;
  281.     clc3   = 1.25 * CK4 * SQR(pInvSQ) * meanMotDeep;
  282.  
  283.     c1SQ   = SQR(c1);
  284.     d2     = 4.0 * sma0Deep * xi * c1SQ;
  285.     clc0   = d2 * xi * c1 / 3.0;
  286.     d3     = (17.0 * sma0Deep + sStar) * clc0;
  287.     d4     = 0.5 * clc0 * sma0Deep * xi * 
  288.              (221.0 * sma0Deep + 31.0 * sStar) * c1;
  289.  
  290.     t2Cof  = 1.5 * c1;
  291.     t3Cof  = d2 + 2.0 * c1SQ;
  292.     t4Cof  = 0.25 * (3.0 * d3 + c1 * (12.0 * d2 + 10.0 * c1SQ));
  293.     t5Cof  = 0.2 * (3.0 * d4 + 12.0 * c1 * d3 + 
  294.              6.0 * SQR(d2) + 15.0 * c1SQ * (2.0 * d2 + c1SQ));
  295.  
  296.     xLCof  = 0.125 * A3OCK2 * sinInc * (3.0 + 5.0 * theta) / (1.0 + theta);
  297.     ayCof  = 0.250 * A3OCK2 * sinInc;
  298.  
  299.     delMeanAn0 = pow(1.0 + eta * cos(meanAn0),3.0);
  300.  
  301.     meanAnDot  = meanMotDeep + 0.5 * clc1 * beta0 * x3th2m1 + 
  302.                  0.0625 * clc2 * beta0 * 
  303.                  (13.0 - 78.0 * theta2 + 137.0 * theta4);
  304.  
  305.     argPerDot  = 0.5 * clc1 * x5th2m1 + 0.0625 * clc2 * 
  306.                  (7.0 - 114.0 * theta2 + 395.0 * theta4) + 
  307.                  clc3 * (3.0 - 36.0 * theta2 + 49.0 * theta4);
  308.  
  309.     raanDot    = (-clc1 - 0.5 * clc2 * x19th2m4 - 2.0 * clc3 * x7th2m3) * theta;
  310.  
  311.     meanAnCof  = -TWOTHIRDS * coef0 * bStar * AE / eccEta;
  312.     argPerCof  = bStar * c3 * cos(argPer0);
  313.     raanCof    = -3.5 * beta02 * clc1 * theta * c1;
  314.  
  315.     return;
  316. }
  317.  
  318. /******************************************************************************/
  319. /*                                                                            */
  320. /* updateGravityAndDrag: updates quantities for secular gravity and           */
  321. /*                       atmospheric drag                                     */
  322. /*                                                                            */
  323. /******************************************************************************/
  324.  
  325. void updateGravityAndDrag()
  326.  
  327. {
  328.     meanAnDF = meanAn0 + meanAnDot * tFP;
  329.     argPerDF = argPer0 + argPerDot * tFP;
  330.     raanDF   = raan0   + raanDot   * tFP;
  331.  
  332.     meanAn   = meanAnDF;
  333.     argPer   = argPerDF;
  334.     raan     = raanDF + raanCof * tSQ;
  335.  
  336.     delA     = 1.0 - c1 * tFP;
  337.     delE     = bStar * c4 * tFP;
  338.     delL     = t2Cof * tSQ;
  339.  
  340.     if (!truncFlag)
  341.     {
  342.         delArgPer  = argPerCof * tFP;
  343.         delMeanAn  = meanAnCof * (pow(1.0 + eta * cos(meanAnDF),3.0) - 
  344.                                   delMeanAn0);
  345.  
  346.         delAM      = delArgPer + delMeanAn;
  347.         meanAn    += delAM;
  348.         argPer    -= delAM;
  349.  
  350.         sinMeanAn0 = sin(meanAn0);
  351.         sinMeanAnP = sin(meanAn);
  352.  
  353.         delA      -= d2 * tSQ + d3 * tCB + d4 * tQD;
  354.         delE      += bStar * c5 * (sinMeanAnP - sinMeanAn0);
  355.         delL      += t3Cof * tCB + t4Cof * tQD + t5Cof * tQN;
  356.     }
  357.  
  358.     smallA  = sma0Deep * SQR(delA);
  359.     smallE  = ecc0 - delE;
  360.  
  361.     xL      = meanAn + argPer + raan + meanMotDeep * delL;
  362.     beta    = sqrt(1.0 - SQR(smallE));
  363.     meanMot = KE / pow(smallA,1.5);
  364.  
  365.     curRaanX       = raan;                                         /* [rad]   */
  366.     curMotionX     = meanMot * MPD / TWOPI;                        /* [rev/d] */
  367.     meanAnomalyX   = meanAn;                                       /* [rad]   */
  368.  
  369.     return;
  370. }
  371.  
  372. /******************************************************************************/
  373. /*                                                                            */
  374. /* updateLongPeridics: updates long periodic quantities                       */
  375. /*                                                                            */
  376. /******************************************************************************/
  377.  
  378. void updateLongPeriodics()
  379.  
  380. {
  381.     upd7  = 1.0 / (smallA * SQR(beta));
  382.  
  383.     axN   = smallE * cos(argPer);
  384.  
  385.     ayNL  = upd7 * ayCof;
  386.     ayN   = smallE * sin(argPer) + ayNL;
  387.  
  388.     xLL   = upd7 * xLCof * axN;
  389.     xLT   = xL + xLL;
  390.  
  391.     curArgPerigeeX = xLT - meanAn - raan;                          /* [rad]   */
  392.     return;
  393. }
  394.  
  395. /******************************************************************************/
  396. /*                                                                            */
  397. /* updateShortPeridics: gets osculating quantities by adding delta terms      */
  398. /*                                                                            */
  399. /******************************************************************************/
  400.  
  401. void updateShortPeriodics()
  402.  
  403. {
  404.     eCosE    = kep4 + kep5;
  405.     eSinE    = kep2 - kep3;
  406.  
  407.     eL2      = SQR(axN) + SQR(ayN);
  408.     upd1     = 1.0 - eL2;
  409.     pL       = smallA * upd1;
  410.     betaL    = sqrt(upd1);
  411.  
  412.     smallR   = smallA * (1.0 - eCosE);
  413.  
  414.     rDot     = KE * sqrt(smallA) / smallR * eSinE;
  415.     rfDot    = KE * sqrt(pL) / smallR;
  416.  
  417.     upd2     = smallA / smallR;
  418.     upd3     = 1.0 / (1.0 + betaL);
  419.  
  420.     cosU     = upd2 * (cosEpAP - axN + ayN * eSinE * upd3);
  421.     sinU     = upd2 * (sinEpAP - ayN - axN * eSinE * upd3);
  422.  
  423.     cos2U    = 2.0 * SQR(cosU) - 1.0;
  424.     sin2U    = 2.0 * sinU * cosU;
  425.  
  426.     smallU   = atan2(sinU,cosU);
  427.     smallU   = reduce(smallU,ZERO,TWOPI);
  428.  
  429.     upd4     = CK2 / SQR(pL);
  430.     upd5     = CK2 / pL;
  431.     upd6     = 1.5 * upd4 * theta;
  432.  
  433.     delR     = 0.5 * upd5 * x1mth2 * cos2U;
  434.     delU     = -0.25 * upd4 * x7th2m1 * sin2U;
  435.  
  436.     delRaan  = upd6 * sin2U;
  437.     delInc   = upd6 * sinInc * cos2U;
  438.  
  439.     delRDot  = -meanMot * upd5 * x1mth2 * sin2U;
  440.     delRFDot = meanMot * upd5 * (x1mth2 * cos2U + 1.5 * x3th2m1);
  441.  
  442.     rk       = smallR * (1.0 - 1.5 * upd4 * betaL * x3th2m1) + delR;
  443.     uk       = smallU + delU;
  444.  
  445.     raanK    = raan  + delRaan;
  446.     incK     = inc0  + delInc;
  447.     rkDot    = rDot  + delRDot;
  448.     rfkDot   = rfDot + delRFDot;
  449.  
  450.     curArgNodeX     = uk;
  451.     curInclinationX = incK;
  452.  
  453.     return;
  454. }
  455.  
  456. /******************************************************************************/
  457. /*                                                                            */
  458. /* keplerX: solves Kepler's equation                                          */
  459. /*                                                                            */
  460. /******************************************************************************/
  461.  
  462. void keplerX()
  463.  
  464. {
  465.     int iterationCount;
  466.  
  467.     capU    = reduce(xLT - raan,ZERO,TWOPI);
  468.     capEpAP = capU;
  469.  
  470.     iterationCount = 0;
  471.  
  472.     do
  473.     {
  474.         kep1     = capEpAP;
  475.         sinEpAP  = sin(kep1);
  476.         cosEpAP  = cos(kep1);
  477.  
  478.         kep2     = axN * sinEpAP;
  479.         kep3     = ayN * cosEpAP;
  480.         kep4     = axN * cosEpAP;
  481.         kep5     = ayN * sinEpAP;
  482.  
  483.         capEpAP  = kep1 + (capU - kep1 + kep2 - kep3) / (1.0 - kep4 - kep5);
  484.         diffAnom = capEpAP - kep1;
  485.  
  486.         iterationCount++;
  487.     }
  488.     while (fabs(diffAnom) >= 1.0e-7);
  489.  
  490.     eccAnom = capEpAP - argPer;
  491.  
  492.     if (fabs(eccAnom) < ONEPPM)
  493.         trueAnomalyX = PI;
  494.     else
  495.         trueAnomalyX = 2.0 * atan(sqrt((1.0 + ecc0) /
  496.                             (1.0 - ecc0)) * tan(eccAnom / 2.0));
  497.  
  498.     trueAnomalyX = reduce(trueAnomalyX,ZERO,TWOPI);
  499.  
  500. /*
  501.     printf("iterations = %d\n",iterationCount);
  502. */
  503.  
  504.     return;
  505. }
  506.  
  507. /******************************************************************************/
  508. /*                                                                            */
  509. /* getOrbitNumberX: gets orbit number and fraction of the orbit               */
  510. /*                  'meanAnomalyX' contains the epoch mean anomaly and the    */
  511. /*                  number of revolutions since the epoch                     */
  512. /*                                                                            */
  513. /******************************************************************************/
  514.  
  515. void getOrbitNumberX()
  516.  
  517. {
  518.     refOrbit   = (double) epochOrbitNum;                             /* [rev] */
  519.     curOrbit   = refOrbit + meanAnomalyX * CRREV;                    /* [rev] */
  520.     orbitNum   = (long) curOrbit + 1L;
  521.     orbitFract = 100.0 * modf(curOrbit,&dummyd);
  522.  
  523.     return;
  524. }
  525.  
  526. /******************************************************************************/
  527. /*                                                                            */
  528. /* getStateVectorX: calculates position and velocity components               */
  529. /*                                                                            */
  530. /******************************************************************************/
  531.  
  532. void getStateVectorX()
  533.  
  534. {
  535.     double mVec[3], nVec[3], uVec[3], vVec[3];
  536.     double cosIK, sinIK, cosRK, sinRK, cosUK, sinUK;
  537.     int    i;
  538.  
  539.     cosUK   = cos(uk);
  540.     sinUK   = sin(uk);
  541.     cosIK   = cos(incK);
  542.     sinIK   = sin(incK);
  543.     cosRK   = cos(raanK);
  544.     sinRK   = sin(raanK);
  545.  
  546.     mVec[0] = -sinRK * cosIK;
  547.     mVec[1] = cosRK * cosIK;
  548.     mVec[2] = sinIK;
  549.  
  550.     nVec[0] = cosRK;
  551.     nVec[1] = sinRK;
  552.     nVec[2] = 0.0;
  553.  
  554.     for (i = 0; i <= 2; i++)
  555.     {
  556.         uVec[i]     = mVec[i] * sinUK + nVec[i] * cosUK;
  557.         vVec[i]     = mVec[i] * cosUK - nVec[i] * sinUK;
  558.  
  559.         satPosS[i]  = rk * uVec[i];
  560.         satVelS[i]  = rkDot * uVec[i] + rfkDot * vVec[i];
  561.  
  562.         satPosS[i] *= EARTHRADIUS;                                  /* [km]   */
  563.         satVelS[i] *= EARTHRADIUS / 60.0;                           /* [km/s] */
  564.  
  565.         satPosGl[i] = satPosS[i];
  566.         satVelGl[i] = satVelS[i];
  567.     }
  568.  
  569.     return;
  570. }
  571.  
  572. /******************************************************************************/
  573. /*                                                                            */
  574. /* calcSGP4: calculates SGP4 propagation model                                */
  575. /*                                                                            */
  576. /* For perigee less than 220 km, the truncFlag is set and the equations for   */
  577. /* 'smallA' and 'xL' are truncated after the 'c1' term, and the terms         */
  578. /*  involving 'c5', 'deltaArgPer' and 'deltaMeanAn' are dropped.              */
  579. /*                                                                            */
  580. /******************************************************************************/
  581.  
  582. void calcSGP4(tmArg)
  583.  
  584. double tmArg;
  585.  
  586. {
  587.     getTimeArgs(tmArg);
  588.  
  589.     if (initNorad)
  590.         calcCommonPerturb();
  591.  
  592.     updateGravityAndDrag();
  593.     updateLongPeriodics();
  594.     keplerX();
  595.     updateShortPeriodics();
  596.     getStateVectorX();
  597.     getOrbitNumberX();
  598.  
  599.     return;
  600. }
  601.  
  602. /******************************************************************************/
  603. /*                                                                            */
  604. /* initDeepSpace:                                                             */
  605. /*                                                                            */
  606. /******************************************************************************/
  607.  
  608. void initDeepSpace()
  609.  
  610. {
  611.     return;
  612. }
  613.  
  614. /******************************************************************************/
  615. /*                                                                            */
  616. /* updateGravity:                                                             */
  617. /*                                                                            */
  618. /******************************************************************************/
  619.  
  620. void updateGravity()
  621.  
  622. {
  623.     return;
  624. }
  625.  
  626. /******************************************************************************/
  627. /*                                                                            */
  628. /* updateDrag:                                                                */
  629. /*                                                                            */
  630. /******************************************************************************/
  631.  
  632. void updateDrag()
  633.  
  634. {
  635.     return;
  636. }
  637.  
  638. /******************************************************************************/
  639. /*                                                                            */
  640. /* deepSpaceSecular:                                                          */
  641. /*                                                                            */
  642. /******************************************************************************/
  643.  
  644. void deepSpaceSecular()
  645.  
  646. {
  647.     return;
  648. }
  649.  
  650. /******************************************************************************/
  651. /*                                                                            */
  652. /* deepSpacePeriodics:                                                        */
  653. /*                                                                            */
  654. /******************************************************************************/
  655.  
  656. void deepSpacePeriodics()
  657.  
  658. {
  659.     return;
  660. }
  661.  
  662. /******************************************************************************/
  663. /*                                                                            */
  664. /* calcSDP4:                                                                  */
  665. /*                                                                            */
  666. /******************************************************************************/
  667.  
  668. void calcSDP4(tmArg)
  669.  
  670. double tmArg;
  671.  
  672. {
  673.     getTimeArgs(tmArg);
  674.  
  675.     if (initNorad)
  676.         calcCommonPerturb();
  677.  
  678.     initDeepSpace();
  679.     updateGravity();
  680.     updateDrag();
  681.     deepSpaceSecular();
  682.     updateLongPeriodics();
  683.     keplerX();
  684.     updateShortPeriodics();
  685.     getStateVectorX();
  686.     getOrbitNumberX();
  687.  
  688.     return;
  689. }
  690.  
  691. /******************************************************************************/
  692. /*                                                                            */
  693. /* getNoradStateVector: calculates state vector following the NORAD SGP4/SDP4 */
  694. /*                      propagation models                                    */
  695. /*                                                                            */
  696. /* uses deep space model if orbital period is longer than 225 minutes         */
  697. /* all time arguments are in minutes                                          */
  698. /*                                                                            */
  699. /* in the SatTrack main code 'getStateVector()' is called before              */
  700. /* 'getNoradStateVector()'                                                    */
  701. /* the former function calculates the 'TLE Mean Model'                        */
  702. /* the latter overwrites the state vector from the former if main() is        */
  703. /* invoked with the '-m' option                                               */
  704. /*                                                                            */
  705. /******************************************************************************/
  706.  
  707. void getNoradStateVector(timeArg)
  708.  
  709. double timeArg;
  710.  
  711. {
  712.     if (modelFlag)                 /* TLE Mean Model will be calculated from  */
  713.         return;                    /* 'getStateVector()'                      */
  714.  
  715.     if (initNorad)
  716.     {
  717.         getSemiMajorAxis();        /* sets 'deepSpaceFlag' if T >= 225 min    */
  718.         checkPerigeeHeight();      /* checks if perigee height < 0            */
  719.     }
  720.  
  721.     if (satCrashFlag)
  722.         return;
  723.  
  724.     if (deepSpaceFlag)
  725.     /*  calcSDP4(timeArg); */      /* overwrites TLE Mean Model state vector  */
  726.         return;                    /* calcSDP4() is still under construction  */
  727.     else
  728.         calcSGP4(timeArg);         /* overwrites TLE Mean Model state vector  */
  729.  
  730.     initNorad = FALSE;
  731.     return;
  732. }
  733.  
  734. /******************************************************************************/
  735. /*                                                                            */
  736. /* End of function block satprop.c                                            */
  737. /*                                                                            */
  738. /******************************************************************************/
  739.